PROGRAM DigiBog_Hydro

!Version: 02/07/2008 (version 2.7)

!-------------------------------------------------------------------------------
! Section 1.0 Program Header.
!-------------------------------------------------------------------------------
!
! DESCRIPTION: 
!   A model to simulate water tables in ombrotrophic bogs using the Boussinesq 
!   equation.  The model uses cm and s.
! 
! CURRENT CODE OWNERS:
!   Paul J. Morris, Andy Baird*, Lisa R. Belyea and Ruth Boogert
!   *School of Geography, 
!   University of Leeds, 
!   Leeds
!   LS2 9JT
!   pmorris@mcmaster.ca
!   a.j.baird@leeds.ac.uk
!   l.belyea@qmul.ac.uk
!   r.boogert@qmul.ac.uk

!  ***** BEGIN LICENSE BLOCK *****
!  Version: MPL 1.1
! 
!  The contents of this file are subject to the Mozilla Public License Version
!  1.1 (the "License"); you may not use this file except in compliance with
!  the License. You may obtain a copy of the License at
!  http://www.mozilla.org/MPL/
! 
!  Software distributed under the License is distributed on an "AS IS" basis,
!  WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
!  for the specific language governing rights and limitations under the
!  License.
! 
!  The Original Code is DigiBog Hydro v2.7.
! 
!  The Initial Developer of the Original Code is
!  Andy Baird & Paul J Morris.
!  Portions created by the Initial Developer are Copyright (C) 2008
!  the Initial Developer. All Rights Reserved.
! 
!  Contributor(s): ----
! 
!  ***** END LICENSE BLOCK *****

! MODIFICATION HISTORY: 
!   Programmer           Date           Modifications
!   ============================================================================
!   Andrew J. Baird      08/04/2005     Original 1.5-D code (v. 1.3)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       24/02/2006     Conversion to 2.5 D and DigiBog Fortran 
!                                       Standards implemented (v. 2.0)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       27/02/2006     Testing, minor correction (v. 2.1)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       05/03/2006     Subroutines 'column_activation' and 
!                                       'steady_state_check' written (v. 2.2)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       19/06/2006     Testing completed (v. 2.3)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       20/03/2007     Code cleaning (v. 2.4)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       24/04/2007     Further cleaning, including removal of
!                                       replicated spatial step reference in
!                                       'move_water' subroutine (v. 2.5)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       09/05/2007     Above-ground storage facilitated in 
!                                       2.5-D version (v. 2.6)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       02/07/2008     Final code cleaning, full annotation
!                                       (v. 2.7)
!   ----------------------------------------------------------------------------
!
! METHOD: 
!   e-mail authors for full documentation.
!    
! INPUT FILES: 
!   All input files are read by the main program, before any subroutines are 
!   called.  These input files must be in the same folder as the project file
!   and all source files.  The input files are:
!
!   hydro_parameters.txt
!     This file contains all the scalar parameters used in the model.
!
!   hydro_no_layers.txt
!     This file contains the number of layers to be used in each vertical, 
!     square-sectioned column.  Sorted in order of x, then y (i.e. for each x 
!     row, values for all y positions are given before moving onto next x).
!
!   hydro_baltitude.txt
!     This file contains the base altitude data for each column, allowing the 
!     user to represent a sloping or uneven substrate if necessary.  Sorted in 
!     order of x, then y.
!
!   hydro_in_wattable
!     This file contains the values for initial water-table heights in each 
!     column.  Sorted in order of x, then y.
!
!   hydro_thickness.txt
!     This file contains the vertical thickness of each peat layer which makes 
!     up the model aquifer.  Sorted in order of x, then y, then z (i.e. one 
!     column at a time, from bottom to top).
!
!   hydro_k.txt
!     This file contains the hydraulic conductivity of each peat layer which 
!     makes up the model aquifer.  Sorted in order of x, then y, then z (i.e. one
!     column at a time, from bottom to top).
!
!   hydro_s_yield.txt
!     This file contains the drainable porosity of each peat layer which makes 
!     up the model aquifer.  Sorted in order of x, then y, then z (i.e. one 
!     column at a time, from bottom to top).
!
! OUTPUT FILES: 
!   Three output files are used, and they are written to at the end of the model
!   run.  
!
!   hydro_wt_output.txt
!     This file records the final water-table position, relative to the 
!     impermeable base, in each column at the end of the model run.
!
!   hydro_x_indices.txt
!     This file records the x position of each column, at the end of the 
!     model run, in order to allow the user to plot the water-table data in an 
!     (x, y, z) triplet.
!
!   hydro_y_indices.txt
!     This file records the y position of each column, at the end of the 
!     model run, in order to allow the user to plot the water-table data in an 
!     (x, y, z) triplet.
!       
! CODE DESCRIPTION: 
!   Language:           Fortran 95 
!   Software Standards: "UK Meteorological Office Fortran 90 Standards".
!   See:
!   http://research.metoffice.gov.uk/research/nwp/numerical/fortran90/
!          f90_standards.html
!
!
! DECLARATIONS 
!
! Modules used:
  USE hydro_procedures

  IMPLICIT NONE 
 
  !Local scalars: 
  INTEGER :: x_extent, & !Model grid extent in x direction                       
             y_extent, & !Model grid extent in y direction   
             z_extent, & !Max number of layers per column                    
             x,        & !Spatial index counter                                   
             y,        & !Spatial index counter                                   
             z           !Spatial index counter
  
  REAL(KIND=q) :: spatial_step,     & !Model horizontal spatial increment
                  timestep,         & !Model temporal increment
                  runtime,          & !Timeout condition
                  elapsed_time,     & !Model clock
                  rainfall,         & !Net rainfall (Precip. - EvapoTrans.)
                  steady_threshold, & !Steady-state criterion per column
                  steady_columns      !Threshold prop'n of cols for steady state

  INTEGER :: alloc_error,   & !Internal array allocation error flag
             active_columns   !Count of active model columns
                  
  CHARACTER(LEN=1) :: param_error      !Internal parameter error flag

  CHARACTER(LEN=9) :: hydro_status     !Indicates model is steady or transient
  
  CHARACTER(LEN=25) :: data_file_name_1,  & !Parameter input file.
                       data_file_name_2,  & !No. layers input file.
                       data_file_name_3,  & !Base altitude input file.
                       data_file_name_4,  & !Initial water table input file.
                       data_file_name_5,  & !Layer thickness input file.
                       data_file_name_6,  & !Layer hydro K input file.
                       data_file_name_7,  & !Layer dr. porosity input file.
                       data_file_name_8,  & !Final WT output file.
                       data_file_name_9,  & !x coordinates output file.
                       data_file_name_10    !y coordinates output file.
      
  !Global/local arrays:

  !activation_status indicates whether column participates in simulation
  CHARACTER(8), ALLOCATABLE, DIMENSION(:,:) :: activation_status
  
  REAL(KIND=q), ALLOCATABLE, DIMENSION(:,:) :: base_altitude, & !Above datum
                                               water_change,  & !
                                               water_table,   & !Above base
                                               wk_mean          !Depth av. K

  !layer_attributes stores layer thickness, hydraulic conductivity and
  !dr. porosity;  transmissivity stores layer elevation and transmissivity;
  !layer_storage stores each layer's water capacity as volume per unit area;
  !i.e. expressed as a depth.                                                   
  REAL(KIND=q), ALLOCATABLE, DIMENSION(:,:,:,:) :: layer_attributes,           &
                                                   transmissivity,             &
                                                   layer_storage 
                                                   
  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: no_layers !Number of layers per column
  

!-------------------------------------------------------------------------------
! Section 2.0 Data Input; Initialisation.
!-------------------------------------------------------------------------------
  
  !Name input data files and output data file
  data_file_name_1  = "hydro_parameters.txt"
  data_file_name_2  = "hydro_no_layers.txt"
  data_file_name_3  = "hydro_baltitude.txt"
  data_file_name_4  = "hydro_in_wattable.txt"
  data_file_name_5  = "hydro_thickness.txt"
  data_file_name_6  = "hydro_k.txt"
  data_file_name_7  = "hydro_s_yield.txt"
  data_file_name_8  = "hydro_wt_output.txt"
  data_file_name_9  = "hydro_x_indices.txt"
  data_file_name_10 = "hydro_y_indices.txt"
  
  !Open files for input of data and file for output of data
  OPEN(UNIT=10,  FILE=data_file_name_1,  STATUS="OLD")
  OPEN(UNIT=20,  FILE=data_file_name_2,  STATUS="OLD")
  OPEN(UNIT=30,  FILE=data_file_name_3,  STATUS="OLD")
  OPEN(UNIT=40,  FILE=data_file_name_4,  STATUS="OLD")
  OPEN(UNIT=50,  FILE=data_file_name_5,  STATUS="OLD")
  OPEN(UNIT=60,  FILE=data_file_name_6,  STATUS="OLD")
  OPEN(UNIT=70,  FILE=data_file_name_7,  STATUS="OLD")
  OPEN(UNIT=80,  FILE=data_file_name_8,  STATUS="REPLACE")
  OPEN(UNIT=90,  FILE=data_file_name_9,  STATUS="REPLACE")
  OPEN(UNIT=100, FILE=data_file_name_10, STATUS="REPLACE")
  
  !Read data from parameter file
  READ (10, *) spatial_step
  READ (10, *) timestep
  READ (10, *) runtime
  READ (10, *) x_extent
  READ (10, *) y_extent
  READ (10, *) z_extent
  READ (10, *) rainfall
  READ (10, *) steady_threshold
  READ (10, *) steady_columns

  !Error check on data
  PRINT *, "The following parameter values have been read from ",              &
                                                                data_file_name_1
  PRINT '(A20, F15.2, A3)', "spatial increment = ", spatial_step,     " cm"
  PRINT '(A20, F15.2, A2)', "time step = ",         timestep,         " s"
  PRINT '(A20, F15.2, A2)', "runtime = ",           runtime,          " s"
  PRINT '(A20, I12, A11)',  "x extent = ",          x_extent,         " columns"
  PRINT '(A20, I12, A11)',  "y extent = ",          y_extent,         " columns"
  PRINT '(A20, I12, A10)',  "z extent = ",          z_extent,         " layers"
  PRINT '(A20, F15.2, A5)', "rainfall = ",          rainfall,         " cm/s"
  PRINT '(A20, F15.2, A3)', "steady threshold = ",  steady_threshold, " cm"
  PRINT '(A20, F15.2)',     "steady columns = ",    steady_columns
  PRINT *, "Are these values correct (Y/N)?"
  READ *, param_error
  SELECT CASE (param_error)
  CASE ("n")
    PRINT *, "Error in data - model terminating"
    STOP
  CASE ("N")
    PRINT *, "Error in data - model terminating"
    STOP
  CASE DEFAULT
    PRINT *, "Data appear to be okay - model run continues"
  END SELECT

!
  !Allocate length to arrays
  ALLOCATE(activation_status(x_extent, y_extent), STAT=alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for activation_status"
    STOP
  END IF
  
  ALLOCATE(no_layers(x_extent, y_extent), STAT=alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for no_layers"
    STOP
  END IF
  
  ALLOCATE(base_altitude(x_extent, y_extent), STAT=alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for base_altitude"
    STOP
  END IF

  ALLOCATE(water_change(x_extent, y_extent), STAT=alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for water_change"
    STOP
  END IF

  ALLOCATE(water_table(x_extent, y_extent), STAT=alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for water_table"
    STOP
  END IF
  
  ALLOCATE(wk_mean(x_extent, y_extent), STAT=alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for wk_mean"
    STOP
  END IF

  ALLOCATE(layer_attributes(x_extent, y_extent, z_extent, 3),                  &
                                                      STAT = alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for layer_attributes"
    STOP
  END IF
  
  ALLOCATE(transmissivity(x_extent, y_extent, z_extent, 2), STAT = alloc_error)
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for transmissivity"
    STOP
  END IF
  
  ALLOCATE(layer_storage(x_extent, y_extent, z_extent, 2), STAT = alloc_error)                                                     
  IF (alloc_error /= 0) THEN
    PRINT *, "Model could not allocate space for layer_storage"
    STOP
  END IF


  !Read data from files to rank-two arrays
  DO x = 1, x_extent
    DO y = 1, y_extent
      READ (20, *) no_layers(x, y)

      !Incorporate above-ground storage layer:
      no_layers (x, y) = no_layers(x, y) + 1
      
      READ (30, *) base_altitude(x, y)
      READ (40, *) water_table(x, y)
    END DO
  END DO

 
  !Initialise layer_attributes:
  layer_attributes = 0.0


  DO x = 1, x_extent
    DO y = 1, y_extent
      
      !Initialise above-ground storage layer properties with the capability to 
      !pond up to 30 cm of surface water:
      layer_attributes(x, y, no_layers(x, y), 1) = 30.0
      layer_attributes(x, y, no_layers(x, y), 2) = 0.0
      layer_attributes(x, y, no_layers(x, y), 3) = 1.0

    END DO
  END DO
  
  !Read data from layer files to rank-four array:
  DO x = 1, x_extent
    DO y = 1, y_extent
      DO z = 1, no_layers(x, y)-1
        
        READ (50, '(3F12.10)') layer_attributes(x, y, z, 1)
        READ (60, '(3F12.10)') layer_attributes(x, y, z, 2)
        READ (70, '(3F12.10)') layer_attributes(x, y, z, 3)
        
      END DO
    END DO
  END DO 
  

  !Initialise remaining arrays:
  transmissivity = 1.0
  layer_storage = 1.0
  wk_mean = 1.0
  water_change = 1.0
  
  !Initialise elapsed time:
  elapsed_time = 0.0



!-------------------------------------------------------------------------------
! Section 3.0 Subroutine Calling; Time Management.
!-------------------------------------------------------------------------------
  !Determine boundary positions and water table levels:
  CALL column_activation(x_extent, y_extent, water_table, activation_status,   &
                         active_columns)


  !Calculate transmissivity profile for each column:
  CALL trans_height(x_extent, y_extent, no_layers, activation_status,          &
                    layer_attributes, transmissivity)
 
  !Calculate amount of water (expressed as a depth) which may be stored in each
  !layer within each column:
  CALL layer_water_depth(x_extent, y_extent, no_layers, activation_status,     &
                         layer_attributes, layer_storage)
  

  !COMMENCE HYDROLOGICAL LOOP
  DO
 
    !(Re-)calculate the depth-averaged K below the water table for each column:
    CALL wat_k_mean(x_extent, y_extent, no_layers, activation_status,          &
                    layer_attributes, transmissivity, wk_mean, water_table)

    !Calculate the amount of water (expressed as a depth) that moves between 
    !each column. The flow law is based on the Boussinesq equation:
    CALL move_water(x_extent, y_extent, activation_status, base_altitude,      &
                    spatial_step, water_table, wk_mean, timestep,              &
                    water_change, rainfall)    

    !Update position of the water table in each column:
    CALL water_table_update(x_extent, y_extent, no_layers, activation_status,  &
                            layer_attributes, layer_storage, water_table,      &
                            water_change)

    !Check for steady-state hydrological behaviour:
    CALL steady_state_check(x_extent, y_extent, activation_status,             &
                            hydro_status, water_change, active_columns,        &
                            steady_threshold, steady_columns)

    !Update elapsed time and check for model termination:
    elapsed_time = elapsed_time + timestep
    IF (elapsed_time > runtime) THEN
      EXIT
    END IF

    !Terminate model run if steady-state has been reached:
    IF (hydro_status == "steady") THEN
      EXIT
    END IF

    !On-screen counter to keep track of time:
    WRITE(*,'(1A4, 20F20.1)') "t = ", elapsed_time

  END DO

!-------------------------------------------------------------------------------
! Section 4.0 Output.
!-------------------------------------------------------------------------------

  !Need to write results to an output file:
  WRITE(80, '(A12)') "WATER TABLES"
  WRITE(80, '(A1)') " "
  DO x = 1, x_extent
    DO y = 1, y_extent
      WRITE (90,  '(1I10)') x
      WRITE (100, '(1I10)') y
      IF(activation_status(x, y) == "off")THEN
        WRITE (80, '(20F20.8)') 0.0
      ELSE
        WRITE (80,'(20F20.8)') water_table(x, y)
      END IF
    END DO
  END DO

    
END PROGRAM DigiBog_Hydro

!-------------------------------------------------------------------------------
! End of Program
!-------------------------------------------------------------------------------
